# Load Packages -----------------------------------------------------------
library(lavaan)
library(haven)

s2 <- read_dta("../Data/study2_raw.dta")

# Merging grid presentation and recoding to pro-trait direction
s2$populism_for_1[which(is.na(s2$populism_for_1))] <- s2$populism_full_1[which(is.na(s2$populism_for_1))]
s2$populism_for_2[which(is.na(s2$populism_for_2))] <- s2$populism_full_2[which(is.na(s2$populism_for_2))]
s2$populism_for_3[which(is.na(s2$populism_for_3))] <- s2$populism_full_3[which(is.na(s2$populism_for_3))]
s2$populism_rev_1[which(is.na(s2$populism_rev_1))] <- s2$populism_full_4[which(is.na(s2$populism_rev_1))]
s2$populism_rev_2[which(is.na(s2$populism_rev_2))] <- s2$populism_full_5[which(is.na(s2$populism_rev_2))]
s2$populism_rev_3[which(is.na(s2$populism_rev_3))] <- s2$populism_full_6[which(is.na(s2$populism_rev_3))]
s2$populism_rev_1 <- s2$populism_rev_1*-1 + 6
s2$populism_rev_2 <- s2$populism_rev_2*-1 + 6
s2$populism_rev_3 <- s2$populism_rev_3*-1 + 6


s2$hs_forward_1[which(is.na(s2$hs_forward_1))] <- s2$hs_full_1[which(is.na(s2$hs_forward_1))]
s2$hs_forward_2[which(is.na(s2$hs_forward_2))] <- s2$hs_full_2[which(is.na(s2$hs_forward_2))]
s2$hs_forward_3[which(is.na(s2$hs_forward_3))] <- s2$hs_full_3[which(is.na(s2$hs_forward_3))]
s2$hs_reversed_1[which(is.na(s2$hs_reversed_1))] <- s2$hs_full_4[which(is.na(s2$hs_reversed_1))]
s2$hs_reversed_2[which(is.na(s2$hs_reversed_2))] <- s2$hs_full_5[which(is.na(s2$hs_reversed_2))]
s2$hs_reversed_3[which(is.na(s2$hs_reversed_3))] <- s2$hs_full_6[which(is.na(s2$hs_reversed_3))]
s2$hs_reversed_1 <- s2$hs_reversed_1*-1 + 6
s2$hs_reversed_2 <- s2$hs_reversed_2*-1 + 6
s2$hs_reversed_3 <- s2$hs_reversed_3*-1 + 6

# Populism-Antidem -------------------------------------------------------------
m <- '
pop =~ populism_for_1 + populism_for_2 + populism_for_3 + populism_rev_1 + populism_rev_2 + populism_rev_3
antidem =~ antidem_1 + antidem_2 + antidem_3 + antidem_4


# error variances
populism_for_1 ~~ populism_for_2 + populism_for_3
populism_for_2 ~~ populism_for_3
'
m_popAntidem <- cfa(m, s2)
summary(m_popAntidem, fit.measures = T, standardized = T)

m_method <- '
# substantive
pop =~ populism_for_1 + populism_for_2 + populism_for_3 + populism_rev_1 + populism_rev_2 + populism_rev_3
antidem =~ antidem_1 + antidem_2 + antidem_3 + antidem_4

# method
method =~ 1*antidem_1 + 1*antidem_2 + 1*antidem_3 + 1*antidem_4 +
1*populism_for_1 + 1*populism_for_2 + 1*populism_for_3 + -1*populism_rev_1 + -1*populism_rev_2 + -1*populism_rev_3

# error variances (content similarities and residual checks)
populism_for_3 ~~ populism_rev_2

# covariances
method ~~ 0*pop + 0*antidem
'
m_popAntidem_method <- cfa(m_method, s2)
summary(m_popAntidem_method, fit.measures = T, standardized = T) 
# adding method factor really improves fit

residuals(m_popAntidem_method, type = "standardized")
modindices(m_popAntidem_method, minimum.value = 10, sort = T)

# store correlations
ests_popAntidem_nomethod <- standardizedSolution(m_popAntidem)
ests_popAntidem_method <- standardizedSolution(m_popAntidem_method)
ests_popAntidem <- rbind(subset(ests_popAntidem_nomethod, op == "~~" & lhs == "pop" & rhs == "antidem", select = c("lhs", "rhs", "est.std", "se", "pvalue")),
                    subset(ests_popAntidem_method, op == "~~" & lhs == "pop" & rhs == "antidem", select = c("lhs", "rhs", "est.std", "se", "pvalue")))
ests_popAntidem$model <- c("No Method", "Method")

# store loadings
load_nomethod <- subset(ests_popAntidem_nomethod, op == "=~", select = c("lhs", "rhs", "est.std", "se"))
load_method <- subset(ests_popAntidem_method, op == "=~" & lhs != "method", select = c("lhs", "rhs", "est.std", "se"))
loadings_popAntidem <- rbind(load_nomethod, load_method)
loadings_popAntidem$model <- c(rep("No Method", nrow(load_nomethod)),
                          rep("Method", nrow(load_nomethod)))
loadings_popAntidem$var <- "popAntidem"

# store error variance
p <- parameterestimates(m_popAntidem_method)
m_var_popAntidem <- p$est[which(p$lhs == "method" & p$rhs == "method")]

# Populism-Conspiracy -------------------------------------------------------------
m <- '
pop =~ populism_for_1 + populism_for_2 + populism_for_3 + populism_rev_1 + populism_rev_2 + populism_rev_3
consp =~ conspre_1 + conspre_2 + conspre_3 + conspre_4


# error variances
populism_for_1 ~~ populism_for_2 + populism_for_3
populism_for_2 ~~ populism_for_3
conspre_1 ~~ conspre_4
'
m_popConspre <- cfa(m, s2)
summary(m_popConspre, fit.measures = T, standardized = T)

residuals(m_popConspre, type = "standardized")
modindices(m_popConspre, minimum.value = 10, sort = T)

m_method <- '
# substantive
pop =~ populism_for_1 + populism_for_2 + populism_for_3 + populism_rev_1 + populism_rev_2 + populism_rev_3
consp =~ conspre_1 + conspre_2 + conspre_3 + conspre_4

# method
method =~ 1*conspre_1 + 1*conspre_2 + 1*conspre_3 + 1*conspre_4 +
1*populism_for_1 + 1*populism_for_2 + 1*populism_for_3 + -1*populism_rev_1 + -1*populism_rev_2 + -1*populism_rev_3

# error variances (content similarities and residual checks)
#populism_for_3 ~~ populism_rev_2
conspre_1 ~~ conspre_4

# covariances
method ~~ 0*pop + 0*consp
'
m_popConspre_method <- cfa(m_method, s2)
summary(m_popConspre_method, fit.measures = T, standardized = T) 
# adding method factor really improves fit

residuals(m_popConspre_method, type = "standardized")
modindices(m_popConspre_method, minimum.value = 10, sort = T)

# store correlations
ests_popConspre_nomethod <- standardizedSolution(m_popConspre)
ests_popConspre_method <- standardizedSolution(m_popConspre_method)
ests_popConspre <- rbind(subset(ests_popConspre_nomethod, op == "~~" & lhs == "pop" & rhs == "consp", select = c("lhs", "rhs", "est.std", "se", "pvalue")),
                         subset(ests_popConspre_method, op == "~~" & lhs == "pop" & rhs == "consp", select = c("lhs", "rhs", "est.std", "se", "pvalue")))
ests_popConspre$model <- c("No Method", "Method")

# store loadings
load_nomethod <- subset(ests_popConspre_nomethod, op == "=~", select = c("lhs", "rhs", "est.std", "se"))
load_method <- subset(ests_popConspre_method, op == "=~" & lhs != "method", select = c("lhs", "rhs", "est.std", "se"))
loadings_popConspre <- rbind(load_nomethod, load_method)
loadings_popConspre$model <- c(rep("No Method", nrow(load_nomethod)),
                               rep("Method", nrow(load_nomethod)))
loadings_popConspre$var <- "popConspre"

# store error variance
p <- parameterestimates(m_popConspre_method)
m_var_popConspre <- p$est[which(p$lhs == "method" & p$rhs == "method")]



# Populism-HS -------------------------------------------------------------
m <- '
pop =~ populism_for_1 + populism_for_2 + populism_for_3 + populism_rev_1 + populism_rev_2 + populism_rev_3
hs =~ hs_forward_1 + hs_forward_2 + hs_forward_3 + hs_reversed_1 + hs_reversed_2 + hs_reversed_3


# error variances
hs_reversed_1 ~~ hs_reversed_2
hs_reversed_2 ~~ hs_reversed_3
populism_for_1 ~~ populism_for_2 + populism_for_3
populism_for_2 ~~ populism_for_3
'
m_popHS <- cfa(m, s2)
summary(m_popHS, fit.measures = T, standardized = T)

m_method <- '
# substantive
pop =~ populism_for_1 + populism_for_2 + populism_for_3 + populism_rev_1 + populism_rev_2 + populism_rev_3
hs =~ hs_forward_1 + hs_forward_2 + hs_forward_3 + hs_reversed_1 + hs_reversed_2 + hs_reversed_3

# method
method =~ 1*hs_forward_1 + 1*hs_forward_2 + 1*hs_forward_3 + -1*hs_reversed_1 + -1*hs_reversed_2 + -1*hs_reversed_3 + 
1*populism_for_1 + 1*populism_for_2 + 1*populism_for_3 + -1*populism_rev_1 + -1*populism_rev_2 + -1*populism_rev_3

# error variances (content similarities and residual checks)
hs_reversed_1 ~~ hs_reversed_2
populism_for_3 ~~ populism_rev_2

# covariances
method ~~ 0*pop + 0*hs
'
m_pophs_method <- cfa(m_method, s2)
summary(m_pophs_method, fit.measures = T, standardized = T) 
# adding method factor really improves fit

residuals(m_pophs_method, type = "standardized")
modindices(m_pophs_method, minimum.value = 10, sort = T)

# store correlations
ests_pophs_nomethod <- standardizedSolution(m_popHS)
ests_pophs_method <- standardizedSolution(m_pophs_method)
ests_popHS <- rbind(subset(ests_pophs_nomethod, op == "~~" & lhs == "pop" & rhs == "hs", select = c("lhs", "rhs", "est.std", "se", "pvalue")),
                     subset(ests_pophs_method, op == "~~" & lhs == "pop" & rhs == "hs", select = c("lhs", "rhs", "est.std", "se", "pvalue")))
ests_popHS$model <- c("No Method", "Method")

# store loadings
load_nomethod <- subset(ests_pophs_nomethod, op == "=~", select = c("lhs", "rhs", "est.std", "se"))
load_method <- subset(ests_pophs_method, op == "=~" & lhs != "method", select = c("lhs", "rhs", "est.std", "se"))
loadings_popHS <- rbind(load_nomethod, load_method)
loadings_popHS$model <- c(rep("No Method", nrow(load_nomethod)),
                           rep("Method", nrow(load_nomethod)))
loadings_popHS$var <- "popHS"

# store error variance
p <- parameterestimates(m_pophs_method)
m_var_popHS <- p$est[which(p$lhs == "method" & p$rhs == "method")]


# Antidem-Conspiracy ------------------------------------------------------
m <- '
antidem =~ antidem_1 + antidem_2 + antidem_3 + antidem_4
consp =~ conspre_1 + conspre_2 + conspre_3 + conspre_4

# error variances. content similarity
conspre_2 ~~ conspre_3
'
m_antidemConspre <- cfa(m, s2)
summary(m_antidemConspre, fit.measures = T, standardized = T)

residuals(m_antidemConspre, type = "standardized")
modindices(m_antidemConspre, minimum.value = 10, sort = T)

m_method <- '
# substantive
antidem =~ antidem_1 + antidem_2 + antidem_3 + antidem_4
consp =~ conspre_1 + conspre_2 + conspre_3 + conspre_4

# method
method =~ 1*conspre_1 + 1*conspre_2 + 1*conspre_3 + 1*conspre_4 +
1*antidem_1 + 1*antidem_2 + 1*antidem_3 + 1*antidem_4

# error variances
conspre_2 ~~ conspre_3

# covariances
method ~~ 0*antidem + 0*consp
'
m_antidemConspre_method <- cfa(m_method, s2)
summary(m_antidemConspre_method, fit.measures = T, standardized = T) 

residuals(m_antidemConspre_method, type = "standardized")
modindices(m_antidemConspre_method, minimum.value = 10, sort = T)

# store correlations
ests_antidemConspre_nomethod <- standardizedSolution(m_antidemConspre)
ests_antidemConspre_method <- standardizedSolution(m_antidemConspre_method)
ests_antidemConspre <- rbind(subset(ests_antidemConspre_nomethod, op == "~~" & lhs == "antidem" & rhs == "consp", select = c("lhs", "rhs", "est.std", "se", "pvalue")),
                         subset(ests_antidemConspre_method, op == "~~" & lhs == "antidem" & rhs == "consp", select = c("lhs", "rhs", "est.std", "se", "pvalue")))
ests_antidemConspre$model <- c("No Method", "Method")

# store loadings
load_nomethod <- subset(ests_antidemConspre_nomethod, op == "=~", select = c("lhs", "rhs", "est.std", "se"))
load_method <- subset(ests_antidemConspre_method, op == "=~" & lhs != "method", select = c("lhs", "rhs", "est.std", "se"))
loadings_antidemConspre <- rbind(load_nomethod, load_method)
loadings_antidemConspre$model <- c(rep("No Method", nrow(load_nomethod)),
                               rep("Method", nrow(load_nomethod)))
loadings_antidemConspre$var <- "antidemConspre"

# store error variance
p <- parameterestimates(m_antidemConspre_method)
m_var_antidemConspre <- p$est[which(p$lhs == "method" & p$rhs == "method")]

# Antidem-HS --------------------------------------------------------------
m <- '
antidem =~ antidem_1 + antidem_2 + antidem_3 + antidem_4
hs =~ hs_forward_1 + hs_forward_2 + hs_forward_3 + hs_reversed_1 + hs_reversed_2 + hs_reversed_3

# correlated errors 
hs_reversed_1 ~~ hs_reversed_2
hs_reversed_2 ~~ hs_reversed_3
'
m_antidemHS <- cfa(m, s2)
summary(m_antidemHS, fit.measures = T, standardized = T)

residuals(m_antidemHS, type = "standardized")
modindices(m_antidemHS, minimum.value = 10, sort = T)

m_method <- '
# substantive
antidem =~ antidem_1 + antidem_2 + antidem_3 + antidem_4
hs =~ hs_forward_1 + hs_forward_2 + hs_forward_3 + hs_reversed_1 + hs_reversed_2 + hs_reversed_3

# method
method =~ 1*hs_forward_1 + 1*hs_forward_2 + 1*hs_forward_3 + -1*hs_reversed_1 + -1*hs_reversed_2 + -1*hs_reversed_3 + 
1*antidem_1 + 1*antidem_2 + 1*antidem_3 + 1*antidem_4

# covariances
method ~~ 0*antidem + 0*hs
'
m_antidemhs_method <- cfa(m_method, s2)
summary(m_antidemhs_method, fit.measures = T, standardized = T)
# fit improves drammatically

residuals(m_antidemhs_method, type = "standardized")
modindices(m_antidemhs_method, minimum.value = 10, sort = T)

# store correlations
ests_antidemhs_nomethod <- standardizedSolution(m_antidemHS)
ests_antidemhs_method <- standardizedSolution(m_antidemhs_method)
ests_antidemHS <- rbind(subset(ests_antidemhs_nomethod, op == "~~" & lhs == "antidem" & rhs == "hs", select = c("lhs", "rhs", "est.std", "se", "pvalue")),
                             subset(ests_antidemhs_method, op == "~~" & lhs == "antidem" & rhs == "hs", select = c("lhs", "rhs", "est.std", "se", "pvalue")))
ests_antidemHS$model <- c("No Method", "Method")

# store loadings
load_nomethod <- subset(ests_antidemhs_nomethod, op == "=~", select = c("lhs", "rhs", "est.std", "se"))
load_method <- subset(ests_antidemhs_method, op == "=~" & lhs != "method", select = c("lhs", "rhs", "est.std", "se"))
loadings_antidemHS <- rbind(load_nomethod, load_method)
loadings_antidemHS$model <- c(rep("No Method", nrow(load_nomethod)),
                                   rep("Method", nrow(load_nomethod)))
loadings_antidemHS$var <- "antidemHS"

# store error variance
p <- parameterestimates(m_antidemhs_method)
m_var_antidemHS <- p$est[which(p$lhs == "method" & p$rhs == "method")]


# Conspiracy-HS -----------------------------------------------------------
m <- '
consp =~ conspre_1 + conspre_2 + conspre_3 + conspre_4
hs =~ hs_forward_1 + hs_forward_2 + hs_forward_3 + hs_reversed_1 + hs_reversed_2 + hs_reversed_3

# error variances. content similarity
conspre_2 ~~ conspre_3
'
m_conspreHS <- cfa(m, s2)
summary(m_conspreHS, fit.measures = T, standardized = T)

residuals(m_conspreHS, type = "standardized")
modindices(m_conspreHS, minimum.value = 10, sort = T)

m_method <- '
# substantive
consp =~ conspre_1 + conspre_2 + conspre_3 + conspre_4
hs =~ hs_forward_1 + hs_forward_2 + hs_forward_3 + hs_reversed_1 + hs_reversed_2 + hs_reversed_3

# method
method =~ 1*hs_forward_1 + 1*hs_forward_2 + 1*hs_forward_3 + -1*hs_reversed_1 + -1*hs_reversed_2 + -1*hs_reversed_3 + 
1*conspre_1 + 1*conspre_2 + 1*conspre_3 + 1*conspre_4

# error variances. content similarity
conspre_2 ~~ conspre_3

# covariances
method ~~ 0*consp + 0*hs
'
m_consprehs_method <- cfa(m_method, s2)
summary(m_consprehs_method, fit.measures = T, standardized = T)
# fit improves drammatically

residuals(m_consprehs_method, type = "standardized")
modindices(m_consprehs_method, minimum.value = 10, sort = T)

# store correlations
ests_consprehs_nomethod <- standardizedSolution(m_conspreHS)
ests_consprehs_method <- standardizedSolution(m_consprehs_method)
ests_conspreHS <- rbind(subset(ests_consprehs_nomethod, op == "~~" & lhs == "consp" & rhs == "hs", select = c("lhs", "rhs", "est.std", "se", "pvalue")),
                        subset(ests_consprehs_method, op == "~~" & lhs == "consp" & rhs == "hs", select = c("lhs", "rhs", "est.std", "se", "pvalue")))
ests_conspreHS$model <- c("No Method", "Method")

# store loadings
load_nomethod <- subset(ests_consprehs_nomethod, op == "=~", select = c("lhs", "rhs", "est.std", "se"))
load_method <- subset(ests_consprehs_method, op == "=~" & lhs != "method", select = c("lhs", "rhs", "est.std", "se"))
loadings_conspreHS <- rbind(load_nomethod, load_method)
loadings_conspreHS$model <- c(rep("No Method", nrow(load_nomethod)),
                              rep("Method", nrow(load_nomethod)))
loadings_conspreHS$var <- "conspreHS"

# store error variance
p <- parameterestimates(m_consprehs_method)
m_var_conspreHS <- p$est[which(p$lhs == "method" & p$rhs == "method")]

# Combine -----------------------------------------------------------------
# correlations
ests <- rbind(ests_antidemConspre, ests_antidemHS, ests_popAntidem, 
              ests_popConspre, ests_popHS, ests_conspreHS)
# method variance
error_vars <- c(m_var_antidemConspre, m_var_antidemHS, m_var_popAntidem, 
                m_var_popConspre, m_var_popHS, m_var_conspreHS)
x <- paste("", round(error_vars, 3))


ests_wide_s2 <- reshape(ests, idvar = c("lhs", "rhs"), timevar = "model", direction = "wide")
ests_wide_s2$study <- "Study 2"

# combine
ests$lhs <- c("Anti-Democratic Attitudes (PW)", "Anti-Democratic Attitudes (PW)",
              "Anti-Democratic Attitudes (PW)", "Anti-Democratic Attitudes (PW)",
              "Populism (B)", "Populism (B)",
              "Populism (B)", "Populism (B)",
              "Populism (B)", "Populism (B)",
              "Conspiracy (PW)", "Conspiracy (PW)")
ests$rhs <- c("Conspiracy (PW)", "Conspiracy (PW)", 
              "Hostile Sexism (B)", "Hostile Sexism (B)",
              "Anti-Democratic Attitudes (PW)", "Anti-Democratic Attitudes (PW)", 
              "Conspiracy (PW)", "Conspiracy (PW)",
              "Hostile Sexism (B)", "Hostile Sexism (B)",
              "Hostile Sexism (B)", "Hostile Sexism (B)")
ests$method_variance <- unlist(strsplit(x, " "))
names(ests) <- c("Var1", "Correlation", "SE", "p_value", "Model", "Var2", "Method Variance")
ests$corr <- paste(ests$Var1, ests$Var2, sep = " with ")
ests_2F <- subset(ests, select = c("corr",  "Correlation", "SE", "p_value", "Model", "Method Variance"))
